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Abstract 

We report on a novel type of instability observed in a noisy oscillator unidirectionally coupled to 
a pacemaker. Using a phase oscillator model, we find that, as the coupling strength is increased, the 
noisy oscillator lags behind the pacemaker more frequently and the phase slip rate increases, which 
may not be observed in averaged phase models such as the Kuramoto model. Investigation of the 
corresponding Fokker-Planck equation enables us to obtain the reentrant transition line between 
the synchronized state and the phase slip state. We verify our theory using the Brusselator model, 
suggesting that this reentrant transition can be found in a wide range of limit cycle oscillators. 
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I. INTRODUCTION 


Synchronization of self-sustained oscillators l|-l3| is crucial in many systems, including 
cardiac cells jd, 5|, circadian clock cells js, 7|, and power grids js-lO]. These systems are 
inevitably subject to various kinds of perturbations such as inherent noise, inhomogeneity 
and environmental changes, and therefore coupling of such oscillators needs to be strong 
enough to overcome the effect of these disturbances and ensure synchronization. 

It is known that strong coupling can be a source of instabilities including oscillation 
death [ll, 12] and chaotic dynamics [l3, However, many of previous studies on coupled 


oscillators in the presence of noise focus on the competition of coupling and noise 
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lid, 


where coupling is expected to suppress noise, leading to fast and stable synchronization. 

In this paper, we present a new synchronization-breaking scenario which can occur in the 
strong coupling regime. Naturally, synchronization does not occur for too weak coupling 
because of the effect of noise. We hnd that, in addition to this trivial desynchronization, 
synchronization is disrupted also for too strong coupling. Such a reentrant transition occurs 
even with very weak noise. We elucidate the condition and mechanism of this reentrant 
synchronization through the analysis of a phase oscillator model. Furthermore, we verify 
that the same reentrant transition occurs in limit-cycle oscillators by using the Brusselator 
model and conhrm the validity of our theory. Our study demonstrates that the reentrant 
transition appears quite generally in coupled noisy oscillators. 


II. COUPLED PHASE OSCILLATOR MODEL UNDER NOISE 

We consider the following phase oscillator that is subject to noise and is influenced by a 
noise-free pacemaker: 

^ = oj + KZ{(t)){h{Vtt) - h{(t))] + ^{t), (1) 

where 0 is the phase, u and 0 are the frequency of the oscillator and the pacemaker, 
respectively, and ^{t) is a Gaussian white noise satisfying = D5{t — t'). Interaction 

is determined by 27r-periodic functions Z{(j)) and h(0). A large class of limit-cycle oscillator 
models can be reduced to this model when the stability of a limit-cycle oscillator is high 
enough compared to noise and coupling strengths j^. Here we adopt the following simple 
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functions: 


Z( 0 ) = sin (0 — a), h{(f)) = —cos (f), ( 2 ) 

with a parameter a. Below we mostly consider the case hi = w. 


A. Averaged model 


Let us first examine the averaged dynamics. When the coupling strength K{> 0) and 
the noise strength D are sufficiently small compared to uj, Eq.([T]) is well approximated by 
an averaged phase model. When hi = cu, the phase difference 'ip = ~ obeys 


ij = KT{ij)+m. 


where the interaction function L is obtained as 


-27r 




27r 


dez{'tp + e){h{d)-h{'tp + e)]. 


The present choice of Z and h yields a Sakaguchi-Kuramoto type interaction function 


(3) 



r('0) = —- {sin('0 — a) + sin a} . (5) 

In the absence of noise, the state -^ = 0 is stable when r'(0) = —|cosa is negative, i.e., 
—7r/2 < a < 7 r/ 2 , which we always consider in the present paper. 

It is convenient to rewrite Eq.(j3]) as 


ip 


dip 




( 6 ) 


where the potential F is given by 


F{iP) ^ - / diP'V{iP') 

Jo 

= — - cos{ip — a) + - (sin a)ip. 


(7) 


li K > D, the phase difference ip tends to stay around the potential minima, and occasionally 
jumps to the two adjacent minima, driven by noise. If a = 0, the right and the left potential 
barriers are the same height, and no net drift appears. On the other hand, if a 7 ^ 0, the 
imbalance of the two barriers causes nonzero drift in the positive direction for a < 0 , or in 
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the negative direction for a > 0. The rate of phase slip is well approximated by Kramers’ 
formula j^: Given a barrier height AF, the rate of overcoming the barrier is proportional 
to exp(—Since the difference between the right and the left barriers is vrsina, the 
net jump rate is proportional to 

/ AFi^K\ f { K'vrsina^ j 

where ATl is the height of the left barrier. That is, in the averaged dynamics, the rate of 
phase slip decreases exponentially as K increases, with its direction determined by the sign 
of a, and no reentrant transition occurs. 

B. Non-averaged model 

When A or D is not sufficiently smaller than cu, the averaging is no longer valid. In 
this case, the situation changes drastically. To see this, we numerically solve Eq. ([1]). The 
frequency oj can be set to unity without loss of generality. Figure d^a) shows three types of 
behavior for different values of K with D = 0.05 and a = 0: For K = 0.01, the system is 
dominated by noise (incoherent region). Synchronization is observed when K is increased 
to K = 1.0. However, as we further increase K xrp to K = 30, synchronization is disturbed 
by a jump of the phase difference by —271 (i.e., a phase slip). Measuring the rate of phase 
slip events against K and D for different values of a, we observe the reentrant transition 
from the synchronized state to the phase slip state [Figd](b)]. The critical value of K for 
this transition decreases as D increases. 

Unidirectionality of the slips is seen in Fig. [I](c). Here, for each K and D, we plot 
(A+ — N_)/ (A_|_ + N_ + e), where A+ and N_ are the number of slips in the positive and the 
negative directions, respectively, and e = 1 is inserted to circumvent zero division. In the 
incoherent region, the slip direction and frequency is determined by the sign of a, reflecting 
the drifting force — ysina appearing in the averaged dynamics. On the other hand, in the 
reentrant region, the slip direction is negative for all a values, implying that the phase slip 
in the non-averaged dynamics is qualitatively different from that in the averaged dynamics. 

The local stability analysis cannot explain the reentrant transition. Linearizing Eq. (lT]) 
around the synchronized state 'ip = 0 with 07 = 0 yields 

'll} =+ C,{t), (9) 
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where \{ut) = Z{ut)h'{ut). The linear stability of the state -0 = 0 is determined by the 
time average of X{u)t), i.e., (a;/27r) X{ut)dt = —(l/2)cosa, which is the same as the 
stablity in the averaged model. The state -0 = 0 is thus linearly stable for — 7 r /2 < a < 7 r/ 2 . 
However, it should be noted that A (cut) can be negative for some range of time even when 
its time average is positive. Hence it seems likely that the system is more easily destabilized 
in the non-averaging model than the averaging one. Nevertheless, the reentrant transition 
is observed even when there is no time interval for the coefficient to be negative. In fact, 
considering the case a = 0, where X{ut) = sin^cut > 0, the synchronized state is never 
destabilized at any time, which suggests that nonlinearity is responsible for the reentrant 
transition. 

To understand the global structure of the system, we again utilize a potential description. 
Equation([T]) can be rewritten as 


'tp 






( 10 ) 


where 


/ dpj'Z+ ut) {hiojt) — hi^p'+ ojt)} 

Jo 

= — - cosi^'ip — a) — - cosi^'ip + 2ut — a) 

1 1 
-I- - cos(2'0 -|- 2ut — a) + 2 

1 1 , 

-I— cos a H— cos(2ci;t — a). 

2 4 


( 11 ) 


( 12 ) 


Note that the potential F{'ip,t) is now time dependent. In general, F{'ip,t) is a 27r-periodic 
function in out. Choosing Z and h as in Eq.(l2]), F is vr-periodic in cut. Figure [2] shows the 
space-time plot of F{'ip,t) for a = 0. It is clearly seen that, in addition to the minimum 
'0 = 0 , which exists in the case of the averaged dynamics as well, there is another minimum 
traveling in the negative direction of tp. One can easily confirm that F{'ip,t) has the traveling 
minimum -0 = —2cut and the two maxima —cut -|- a and —cut + a + ii. It is expected that, in 
the non-averaged dynamics, phase slip occurs along this traveling minimum. 

The potential is not tilted if a = 0, and tilted if a 7 ^ 0. Since the former is easier to 
analyze, below we first present our analysis for a = 0 , and then extend the analysis to the 
case a ^ 0. 
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FIG. 1. (Color online) (a) Time evolution of the phase difference '0 for three different values of 
K, with a = 0 and D = 0.05. (b) Phase slip rate against the coupling strength K and the noise 
intensity D ioY a = 0, 0.3, and —0.3. Theoretical lines are given by Eas. lfTSl) and (fT^ for K >1. (c) 
Unidirectionality of slip against K and D. Red and green indicate positive and negative directions 
of '0, respectively. 

1. Non-tilted potential (a = Q) 


The mechanism and condition of the phase slip can be understood through the analysis 
of the following Fokker-Planck equation; 




dt 


dF 


dif] \d'ip 


P\ + 


Dd'^P 

2 


(13) 


where P{'ijj,t) is the probability distribution function for the phase difference ip. By numer¬ 
ically solving Eq. flldh with Eq. flT^ for a = 0, we obtain time evolution of P{'ip,t). Figure 
n a) shows P{'ip,t) together with the trajectories of the extrema of F{'ip,t). Note that, for 
a = 0, three of them cross each other at the same time. When the coupling is sufficiently 
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FIG. 2. (Color online) Space-time plot of the time dependent potential (Eq.(([T2]))) for 

a = 0. 

large, P splits into two components at some moment, one localized at ^ = 0 corresponding 
to the synchronized state, and the other traveling along '0 = —2ut, which corresponds to 
the phase slip state. This traveling component appears only for sufficiently strong coupling, 
although the distribution around -^ = 0 is sharper for stronger coupling [Fig|3](b),(c)]. 

The scenario of how the traveling component emerges is as follows. Let us focus on a 
short time interval around t = 0. When t < 0, P is localized at "0 = 0. At t = 0, the 
three extrema of F cross each other [Fig. [3](a)]. Around this time, the curvature of the 
potential at = 0 almost vanishes, and hence the diffusion dominates the dynamics. Then, 
for t > 0, a potential maximum located at ip = —ujt gradually grows, and, at some time 
tc the drift force caused by the potential becomes comparable to the effect of diffusion. At 
this moment, the part of P{ip,t) located beyond the maximum (i.e., ip < —ojt) is separated 
from the component around ip = 0 and thus conveyed with the potential minimum located 
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at -0 = —2ujt. This process repeats itself with period vr. 

We can roughly estimate the above tc and the transition line through a dimensional 
analysis. First, tc is determined by balancing the drift and the diffusion terms in Eq. fll3p 
evaluated at the potential maximum '0 = —cut: 

KF{-utc,tc) D. (14) 

Because diffusion is dominant for 0 < t < tc, the width of P grows roughly as '/Dt within 
this duration. If this width is comparable to distance to the potential maximum at t = tc, 

i.e., 

^/D^c ~ cute, (15) 

a substantial part of P will be conveyed. From Eqs. flTTll and flTHll . and noting that the 
potential height for small t is F(—cut,t) ~ (cut)"^, we obtain the following scaling relation: 

Il~cu5jF“T (16) 

By substituting Eq. flTB]) into Eq. flTSD . we obtain the relationship between K and tc'. 

tc^u!~^K~^. (17) 

The reason why stronger coupling induces more phase slips is now clear: Since larger coupling 
strength K implies smaller tc, the condition y/Dp > cute is easier to satisfy. 


2. Tilted potential (a ^ 0) 


In the case of small but nonzero |q:|, where the three extrema cross each other at different 
timings, a similar argument can still be made to obtain the transition lines: If a > 0, ^ = 0 
becomes unstable when p = —2cut crosses 0 at t = 0. On the other hand, if a < 0, -0 = 0 
loses its stability at t = —cu“^|a|, when p = —cut + a crosses '0 = 0. In both cases, the 
potential barrier is located at 0 = —cut + a. These situations are schematically shown in 
Fig. m Hence, instead of Eq. flT^ . we have different balance equations Xy/Dp = ojtc — a 
{a > 0) and Xy/D{p — cu ^a) = ujtc — a (a < 0), where a parameter A has been introduced. 
The transition line is thus obtained by eliminating tc from the following equations: 


C 


KP{-tc + a,tc) = D, 

X‘^oj~^D + a (a < 0) 
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(18) 


cut. 


(a > 0), 


(19) 









Substituting a = 0 obviously reproduces Eqs. flT^ and flT5|) . Theoretical lines given by 
Eqs. fllSp and flT^ with A = 3.0 agree well with numerical data for a = 0 and a = 0.3 in 
Fig. [T]^b), although for a = —0.3 it slightly deviates from the numerical data particularly for 
large K. This result suggests that the sign of a is critial for stable synchronization; phase 
slip hardly occurs if a < 0 even for large K. 

3. Frequency mismatch 

We can further extend our scaling argument for a = 0 to the case VL ^ uj. Since the 
frequency difference Acu = a; — serves as a drifting force, Eq. lfTSj) is now modihed as 

X\/ Dtc — (J Acute = 1^4, (20) 

where a > 0 is another parameter. Combining Eq. fl20l) with Eq. flTTll . we obtain 

D ~ CU3 + (T 

Figure [5](a) shows the phase slip rates for different values of Acu. Rescaling according to 
Eq. fl2l]) results in the collapse of data points [Fig|5](b)], which strongly supports the validity 
of our theory. 

C. Generality of reentrant transition 

We have studied the non-averaged phase model ([I]) using a specific form of Z and h as 
in Eq.(l2]). Now let us consider reentrant transition in more general settings. For simplicity, 
we again focus on the case hi = cu. As we have observed above, the extrema of the potential 
F('0,t), which are given by zeroes of + ut) and Hit'll),t) = h{ujt) — + cut) for 

0 < "0 A 27r, play important roles on phase slip. Obviously H has a zero at ^ = 0, 
which corresponds to the synchronized state. Now we assume that the functions Z and 
h are unimodal and that Z has two zeroes, which are reasonable in many practical cases. 
Then it follows that H has another zero ^ = /S(t), which we call the moving TT-branch, 
corresponding to phase slip. In fact, from the above assumption it is easy to show that 
Pit) is a monotonically decreasing function of t that satishes < —cu, crossing ip = 2mi 
(u = 0,1, 2, • • •) twice within one oscillation period, reflecting the fact that phase slips occur 


cu / 


( 21 ) 
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in the negative direction. In addition, two zeros of the response function Z^ijj + ut), denoted 
by ai and a 2 , give the other potential extrema moving linearly along -0 = —ujt + ai^ 2 , which 
we call Z-branches. Note that, in the case of Eq.([2]), the Z-branches are 0 = —ut + a and 
—ut + a + IT, and the moving H branch is 0 = —2ut. 

The synchronized state 0 = 0 is destabilized by the crossing of the Z-branches and the 
moving //-branches, and phase slip may occur, depending on in which order these branches 
cross 0 = 0: if the branch crossing is like FigJU^a), that is, if 0 = 0 is destabilized by the 
moving H branch Erst, and then re-stabilized by a Z-branch, phase slip is more likely to 
occur; if the crossing is like FigJU^b), phase slip is less likely. 


III. REENTRANT TRANSITION IN COUPLED LIMIT CYCLE OSCILLATORS 

Let us conhrm that the reentrant transition occurs in limit cycle oscillators. We demon¬ 
strate it by using the Brusselator model: 

Ui = A - {B + l)ui + u^Vi + K'f’iuj - Ui) + iiit), (22) 

Vi = Bui - u^Vi + //fVi - u) + (23) 

where (0 j) = (1,2) or (2,1), Uij and Vij are the state variables, and ^i(t) and rjiit) are the 
Gaussian white noise satisfying {^iit)^i{t')) = {f]i{t)r]i{t')) = DiS(t — t'). Here we consider 
the case where oscillator 2 acts as a noise-free pacemaker: D 2 = 0, = 0. On 

the other hand, oscillator 1 is under noise, Di = D, and is influenced by oscillator 2 via 
M-coupling (/l|“^ = //, = 0) or u-coupling = 0, = K). 

In the case of u-coupling, too strong coupling results in occasional oscillation failures, 
as shown in Figj6](a). If we measure the phase of the state (m,u) in a standard way j^, 
oscillation failure can be expressed as phase slip. The frequency of phase slips increases as 
K increases for fixed D, indicating a reentrant transition from synchronization to phase slip 
[Figj6](b)]. On the other hand, u-coupling does not induce the reentrant transition [Figj6](c)]. 

The difference between the two ways of coupling can be understood by applying the 
analysis in Sec. Ill Cl to the phase model of the Brusselator. We numerically solve Eq.([I]) 
with Z and h calculated from the numerical phase reduction of Eqs. fl2^ and fl23l) . The 
trajectory of the phase difference 0(f) and the potential extrema are shown in Fig|71 Here 
we can clearly see the two types of branch crossing discussed in Sec. Ill Cl In the case of 
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u-coupling [FigJ3^a)], it is observed that occasionally leaves the state 'ip = 0 and travels 
along the moving if-branch. This occurs when the state ip = 0 becomes temporally unstable 
due to the crossing by the moving ii-branch. Its stability recovers when the state = 0 is 
crossed by one of the Z-branches. This branch crossing is the type a > 0 [FigJl](a)]. On the 
other hand, in the case of w-coupling, phase slip is not observed [FigJTl^b)], where branch 
crossing is of the type a < 0 [FigJU^b)]. These results correspond to the presence and the 
absence of the reentrant transition in Figj6](b) and (c), respectively. 

It is remarkable that the reentrant transition in the Brusselator model is well understood 
from the analysis of its corresponding phase model for such a strong coupling case. Our 
results indicate that non-averaged phase models are useful for the understanding of strongly 
coupled oscillators. 


IV. CONCLUDING REMARKS 


Our study has uncovered that, in addition to a lower limit in coupling strength, there 
is generally an upper limit over which synchronization is disrupted by phase slips when 
an oscillator is subject to noise. Therefore, when strong coupling is required, the way of 
coupling should carefully be constructed. 

In this study we have only considered unidirectional coupling. In the case of mutual 
coupling, we have confirmed in a preliminary numerical study that mutual coupling also 
exhibits reentrant transition, although oscillation death is more commonly observed. 

A similar eentrant transition is known to occur in a certain class of chaotic oscillators 



Our study would further motivate such studies. 
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FIG. 3. (Color online) (a) Space-time plot of the probability distribution function P{ijj,t) obtained 
by solving Ea. (ll3p with K = 30, D = 0.1, and a = 0. Only the region > 0.01 are shown 

(red). The minima (solid lines) and the maxima (dotted lines) of the potential F('0,t) are also 
shown, (b, c) Snapshots of P(V’,t) for K = 1 (dotted) and K = 2>Q (solid), with D = 0.1 and 
a = 0. Shaded regions represent (b) t = 0, (c) t = 1.25. Probability current corresponding 

to the phase slip is indicated by an arrow in (c). 
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FIG. 4. Schematic of potential extrema crossing for nonzero a: (a)a > 0, (b)a < 0. Solid (dotted) 
lines correspond to minima (maxima) of the potential 


(a) 



(b) 



FIG. 5. (Color online) (a) Phase slip rates against K for different w = Aca + fl, with D = 0.3, 
a = 0 and = 1. (b) Same data plot for larger K, scaled by (1 + crAw/o;)®, with the fitting 
parameter cj=1.35. 
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FIG. 6. (Color online) Two nnidirectionally coupled Brusselators: Parameters are A = 1.6 and 
B = 5.0. (a) Time series of ui{t) (red) and U 2 {t) (bine) in the case of 'c-coupling, with K = 0.5 
and D = 0.01. Oscillation failures are observed at t ~ 220 and t ~ 280, indicated by arrows, (b, 
c) Frequency of phase slip against D and K: (b)n-coupling, (c)n-coupling. 
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FIG. 7. (Color online) Trajectories (gray) of the phase difference ^p{t) obtained by solving Eq.([T]), 
with Z and h calculated from the numerical phase reduction of Eas. (|22p and (I23p . with D = 0.005, 
K = 1.0. The extrema of the numerically constructed potential F{'tp,t) are shown together: ip = 0 
(green), Z-branches (red), and a moving 77-branch (blue). Solid and dotted lines correspond to 
minima and maxima of the potential, respectively. The trajectory is plotted modulo 27r in cut. 
(a)u-coupling, (b)u-coupling. 
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